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Na-montmorillonite hydrates in presence of ethane molecules are studied by means of hybrid 
Monte Carlo simulations in the NP ZZ T and fJ,P zz T ensembles. The NP ZZ T ensemble allows us to 
study the interlaminar distance as a function of water and ethane content. These data show clear 
plateaus for lower ethane contents and mainly for water contents consistent with the formation of 
a single water layer. In addition, from this ensemble the structure for some of these interlaminar 
compositions were analyzed. For systems containing few ethane molecules and water enough to 
complete a single layer, it was observed that ethane mainly situates close to the interlayer midplane 
and adopts a nearly parallel arrangement to the clay surface. On the other hand, the [iP zz T 
ensemble allows us to determine the interlaminar distance and water-ethane content for any specific 
reservoir. Here, some important findings are the following: the partial exchange of water by ethane 
molecules that enhances for decreasing the water vapor pressure; the obtention of a practically 
constant interlaminar space distance as a function of the water vapor pressure; the conservation of 
ion solvation shells; the enhancement of the water-ethane exchange for burial conditions; and finally, 
the incapability for a dehydrated clay mineral to swell in a dry and rich ethane atmosphere. 

PACS numbers: 



I. INTRODUCTION 

Beyond the role that clay minerals play in several en- 
vironmental and industrial processes, one major concern 
is how their presence affect petroleum migration and oil- 
recovery |g. In these processes, petroleum, natural 
gas, and their associated brines interact with clay min- 
erals leading to retention, release, and selective adsorp- 
tion of water, ions, and certain organic species, affecting 
the formation permeability. In general, the amount and 
kind of the adsorbed species depend on their chemical 
potential, which in turn establishes the swelling of a spe- 
cific clay, for a given pressure and temperature. Hence, 
changes in the local reservoir composition or conditions, 
due to, for example, an oil-recovery process, can cause 
a destabilization of the formation clay minerals. This 
destabilization, together with a high percentage of clay 
mineral in the basin, may lead to a formation damage, 
i. e. diminution of the permeability. Therefore, knowl- 
edge of the clay behavior in presence of petroleum and 
natural gas compounds and its associated brines at reser- 
voir conditions is mandatory for a better understanding 
of such processes. 

One possible way to study these systems is by means 
of computer simulations. In the last few years, this 
tool has been widely used for studying the swelling of 
montmorillonite hydrates for different interlayer cations 
SSllSSSIIimm There are, however, 
few theoretical works devoted to study the montmoril- 



* godriozo@imp.mx 
1 aguilarf@imp.mx 
-! jllemus@imp.mx 



lonite behavior under basin conditions 0, Q, , and 
there are also few theoretical papers dealing with the 
interaction of montmorillonite hydrates with petroleum 
species Hence, questions such as do any 

oil or natural gas compounds enter into the interlami- 
nar space at reservoir conditions?, how many and what 
kind of molecules enter?, do they expel water molecules 
in the process?, and how does this affect the interlaminar 
space?, still do not have clear answers. 

For tackling some of these questions we performed hy- 
brid Monte Carlo (HMC) simulations to study the behav- 
ior of Na-montmorillonite hydrates at equilibrium with 
different ethane-water reservoirs. Ethane was chosen 
simply because it is an abundant natural gas and live-oil 
compound, and it is expected to enter into the interlam- 
inar space. In addition, the Na + ion was chosen since it 
is the most abundant cation found in cation-clay-water 
systems. Simulations were carried out in the NP ZZ T and 
fiP zz T ensembles. This last one has the advantage of al- 
lowing a direct measurement of the amount of molecules 
and the interlaminar distance at equilibrium under the 
established conditions. The employed water model is 
based on the TIP4P- water |2(j and the Na + -water-clay 
interactions are given by Boek et al. l2ll. The ethane 
parameters were taken from Nath et al. |22| . Two condi- 
tion sets were studied. These are ground level, P= 1 atm 
and T= 298K, and P= 600 atm and T= 394K, which 
corresponds to a burial depth of 4 km (assuming average 
gradients of 30 K/km and 150 atm/km). 

The paper is organized as follows. In Sec.[n] we briefly 
describe the models and the methodology employed for 
performing the simulations. The results are shown in 
Sec. IIIII Finally, Sec. IIVI summarizes the main results 
and extracts some conclusions. 
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II. METHODOLOGY 



A. The model 



The Wyoming type montmorillonite cell given by Skip- 
per et al. [23| was replicated in such a way to form a 4 x 2 
layer. This layer has lateral dimensions of L x = 21.12 A 
and L y = 18.28 A, and thickness of 6.56 A. To obtain 
a typical Wyoming type montmorillonite, isomorphous 
substitutions of trivalent Al atoms of the octahedral sites 
by divalent Mg atoms, and tetravalent Si by trivalent 
Al atoms were made. In this way, Clay I of the work 
by Chavez-Paez et al. |24j| was obtained, with unit cell 
formula Na . 75 nH2O(Si 7 . 75 Al .25)(Al3.5Mg . 5 )O 20 (OH)4. 
Two layers were considered in the simulation box to avoid 
system size effects 1241. The rigid TIP4P model was used 
for water molecules j20| and a flexible model was consid- 
ered for ethane molecules, just as described in the work 
by Nath et al. 22] . For the initial configuration, water 
and ethane molecules were randomly placed in the in- 
tcrlaminar spaces, and sodium ions were distributed in 
the interlayer midplanes. Note that six sodium ions per 
interlaminar space are needed to keep the system elec- 
troneutral. Periodical boundary conditions were applied 
in the three space directions. 

The site-site intermolecular interactions are given by a 
Coulombian contribution plus a Lennard-Jones potential, 
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where subindices i and j are molecular indexes, and a 
and b run over all sites of each molecule. Here, q a and qb 
are the charges of the corresponding sites, e a b and a a b are 
site-to-site specific Lennard-Jones parameters, and r a b is 
the inter-site distance. The Lennard-Jones parameters 
for single sites are given in Table Q] The bond stretching 
potential for ethane molecules is U(r)—K r (r — b eq ) 2 /2, 
with A' r =191.765 kcal/(mol A) and & e9 = 1.54 A The 
site to site Lennard-Jones parameters are given by the 
Lorentz-Berthelot rules 
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We should point out that the oxygen-oxygen parame- 
ters' values were chosen to be exactly those of the TIP4P 
water model. These values were slightly changed in ref- 
erences HH 113 • In addition, the sodium ion param- 
eters were obtained by a fitting procedure that leads to 
a relatively good agreement between the pair energies 
obtained from equations (7) and (8) of Bounds's work 
[2^1 and equation These are, in fact, also similar to 
those reported by Boek et al. |2l| . The fitting procedure 
was performed by employing the Levenberg-Marquardt 



TABLE I: Lennard-Jones parameters for ethane-water-clay- 
Na + interactions. 



Sites 



: (kcal/mol) 



a (A) 



o 

H 

Na 
Si 
Al 
Mg 
CH 3 



0.155 
0.000 
1.570 
3.153 
3.153 
3.153 
0.1990 



3.154 
0.000 
1.740 
1.840 
1.840 
1.840 
3.825 



algorithm and considering several positions and orienta- 
tions of the water- ion pair. To check the accuracy of 
the fitted expression, a HMC simulation of 216 water 
molecules and a Na + ion (and a CI - ion to keep the sys- 
tem electroneutral) was performed. The obtained Na + - 
oxygen and Na + -hydrogen radial distribution functions, 
g(r), and coordination numbers, n(r), were very similar 
to those reported by Bounds 26] . Finally, parameters for 
Si were taken from Marry et al. j^- Parameters for Al 
and Mg were assumed to be equal to those of Si. 

The Ewald summation formalism was implemented for 
handling the electrostatic interactions. For this purpose, 
we set the convergence factor to 5.6/ L m i n , being L m i n 
the minimum simulation box side. Five reciprocal lattice 
vectors for the directions along the shortest sides and 
six for the direction along the largest side were set [28| . 
For the Lennard-Jones contribution a spherical cutoff of 
L min /2 was considered. In addition, this contribution 
was corrected by the standard methods for homogeneous 
fluids IH. 



B. Simulations 

A hybrid Monte Carlo scheme was employed for the 
simulations. This technique has the advantage of allow- 
ing global moves while keeping a high average acceptance 
probability. A global move in configuration space con- 
sists in assigning velocities and integrating the system 
over phase space using some discretization scheme |3(j . 
Velocities are randomly assigned from a Gaussian distri- 
bution in correspondence with the imposed temperature 
and in such a way that total momentum equals zero for 
both interlaminar spaces. The discretization scheme we 
employed is a reversible multiple time scale algorithm 
based on the Trotter expansion of the Liouville propaga- 
tor, which has the advantage of allowing to split forces 
into short and long range ones [3l]]. Hence, a short time 
step is employed to generate the motion using the rapid 
varying short range forces, and a long time step is used 
only to correct the velocities by using the slow varying 
long range forces. The short range forces were chosen as 
the Lennard-Jones contribution plus the real part of the 
electrostatic potential, and the long range forces were set 



3 



as the reciprocal space contribution of the electrostatic 
potential. This allows reducing the expensive calls to 
the reciprocal space contribution. For diminishing time 
correlations, a new configuration is obtained after ten in- 
tegration steps. Once the new configuration is generated, 
it is accepted with a probability 



P = min{l,exp(-/3AH)} 



(4) 



where ATL is the difference between the Hamiltonians 
associated to the new and previous configuration and (3 
is the inverse of the thermal energy. Note that increasing 
the time step will result in a larger deviation from total 
energy conservation and so, in a lower average acceptance 
probability. We fixed the long time step to eight times 
the short time step, and we set the short time step to 
obtain an average acceptance probability of 0.7 [sTj] ■ For 
more details about the hybrid Monte Carlo algorithm see 
Mehlig et al. H3. 

For sampling in the NP ZZ T ensemble, after a trial 
change of particles' positions, a box change is attempted 
in such a way that the stress normal to the surface of 
the clays, P zz , is kept constant. For this purpose, box 
fluctuations are allowed only in the z-direction and the 
probability for accepting a new box configuration is given 
by 

P=min{l,exp[-/3(AW + P ZZ AV- Np-^VjVo))]} 

(5) 

Here, AU is the change in the potential energy, AV is the 
volume change, N is the total number of molecules, and 
V n and V are the new and old box volumes, respectively 

For sampling in an open ensemble, the possibility of 
insertions and deletions of water and ethane molecules 
has to be considered. Water insertions and deletions 



were performed by Rosenbluth sampling as explained 
elsewhere The only difference is that we account 

for the real part of the electrostatic contribution for the 
hydrogen trial conformations (including the non massive 
TIP4P site contribution). This is a necessity since the 
TIP4P hydrogen sites do not have a dispersion-repulsion 
contribution, in contrast to the MCY water model they 
employed. Hence, the real part of the electrostatic con- 
tribution appears directly in the Rosenbluth factors and 
must not be included in the exponential part of equations 
(5) and (6) of the work by Hensen et al. j^. Similarly, 
insertions and deletions of ethane molecules were also 
performed by the Rosenbluth method. Here, the slight 
difference is that ethane molecules are flexible. In this 
case, the method is clearly explained elsewhere [33^. 

Commonly, the open ensemble employed in simulations 
is the Grand canonical, i. e. the fiVT ensemble. This was 
successfully employed to determine the stable hydration 
states of the water-ion-clay systems, i. e. the number 
of water molecules and the interlaminar space at equi- 
librium with a given reservoir 0, Q 0, H3> EH • For 
this purpose, it is necessary to measure the pressure as 
a function of the interlaminar distance (volume) to de- 
termine the states that produce the reservoir's pressure. 
The problem is that this implies a large number of sim- 
ulation runs. Another possibility would be sampling in 
the \xP zz T ensemble. In this way, a single simulation run 
is employed for obtaining the system at equilibrium with 
the reservoir's conditions. 

For finding an algorithm capable for sampling such en- 
semble, it is necessary to pay attention on the probability 
density of finding the system in a particular configura- 
tion. This is deduced from the partition function, which 
in this case reads as 
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where A is the thermal de Broglie wavelength, subindexes 
1 and 2 refer to components 1 and 2 (water and ethane), 
and s N are scaled coordinates 33] . Hence, the corre- 
sponding probability density of finding the system in a 
particular configuration is given by 



V N cxp{-/3[U(s N )-^ 1 N 1 -n 2 N 2 + P zz V}} 
A MiA12 p^t oc ———^ 

(7) 

and so, the algorithm must sample this distribution. For 
this purpose, it is easy to show that particle movements, 
insertions and deletions, and box changes must be done 
as in typical NVT, fiVT and NP ZZ T sampling H3- In 



particular, after trying a change of particles' positions, 
we perform several attempts of inserting-deleting water 
and ethane molecules. This is done by randomly calling 
the four possible trials in such a way that all calls are 
equally probable. This is important in order to guaranty 
the establishment of detail balance. Since accepting in- 
sertions or deletions are rare, we repeat this step 10 times 
or until accepting any insertion or deletion. In case of re- 
fusing the 10 insertion-deletion trials, we performed a box 
trial move. In this way, the system rapidly evolves to an 
equilibrium state that, in general, depends on the initial 
conditions. 
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FIG. 1: Interlaminar distance as a function of the number 
of water molecules per interlaminar space. Symbols □, o, a, 
v, o and < correspond to 0, 6, 9, 12, 15 and 18 fixed ethane 
molecules, respectively. 



FIG. 2: Interlaminar distance as a function of the number 
of ethane molecules per interlaminar space. Symbols □, o, 
a, v and o correspond to 0, 16, 32, 48 and 62 fixed water 
molecules, respectively. 



III. RESULTS 

In order to organize the results this section is split into 
two parts. These are Sampling in the NP ZZ T ensemble 
and Sampling in the [iP zz T ensemble. Each part presents 
the results for the conditions of ground level and for 4 km 
of burial depth. These are T= 298 K, P= 1 atm and T= 
394 K, P= 600 atm, respectively. 



A. Sampling in the NP ZZ T ensemble 

This ensemble allows measuring the interlaminar dis- 
tance as a function of the imposed number of water and 
ethane molecules for a given temperature and pressure. 
This means that water and ethane content may or may 
not correspond to real systems. 

Fig- D shows the interlaminar distance as a function of 
the number of water molecules per clay sheet. Here, the 
ethane content is fixed and the established conditions 
correspond to ground level. Naturally, it is observed 
that as far as the ethane content increases, the curves 
reach larger interlaminar distances. Hence, the bottom 
curve corresponds to zero ethane molecules and the up- 
permost curve to 18 ethane molecules per interlaminar 
space. It should be pointed out that this lowest curve is 
in good agreement with others reported elsewhere [24| . 
suggesting that the differences in models and methods 
are not very important. This curve starts from an in- 
terlaminar distance of 10.22 A, where ions are the only 
specie. This dehydrated state is in good agreement with 
those reported for the TIP4P water model ^HMUl and 
slightly larger than the reported experimental values (9.8 
- 10.0 A) |35| . The curve grows with increasing water con- 



tent reaching a plateau for 20-48 water molecules. This 
plateau yields interlaminar spaces in the range 12.0-12.7 
A, in good agreement with experimental data. For larger 
water contents, a second plateau is developed, which cor- 
responds to a double water layer. 

When six ethane molecules are added to the interlam- 
inar space, a large plateau is developed in the water 
content range of 0-27, which corresponds to interlami- 
nar spaces ranging in 12.2-12.5 A. This means that as 
early as few ethane molecules enter into the interlaminar 
space, they force it to yield relatively large interlaminar 
distances and so, additional water molecules easily en- 
ter without changing it. Furthermore, the fact that only 
27 + 6 total molecules saturate the first layer implies that 
ethane molecules occupy approximately the triple effec- 
tive volume than water. This is verified by the fact that 
for 9 and 12 ethane molecules the plateaus shorten, dis- 
appearing for 15 ethane molecules per interlaminar space. 
Note that in each case the plateaus are yielded for 0-18 
and 0-6 water molecules, respectively. We should also 
pointed out that in all cases the plateaus produce inter- 
laminar spaces in the range of 12.2-12.5 A. Hence, if 
ethane goes into the interlaminar space, it is expected to 
displace water in the process without affecting the inter- 
laminar distance. In addition, these facts make us think 
that for a single layer and in the most favorable situa- 
tion, 12 ethane molecules may enter into the interlami- 
nar space, displacing approximately 36 water molecules 
in the process. 

Similar conclusions can be drawn by studying the in- 
terlaminar distance as a function of the number of ethane 
molecules keeping fixed the number of water molecules. 
This was done and it is presented in Fig. |3 As can be 
seen, the plateaus for 0,16 and 32 water molecules corre- 
spond to the ranges 6-16, 3-9 and 0-3 ethane molecules, 
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FIG. 3: a) Density profiles of oxygen, hydrogen, methyl, and 
sodium sites. They were obtained by fixing 27 water and 
6 ethane molecules, b) The corresponding second Legendre 
polynomial order parameters for water and ethane (see text). 

respectively. In addition, the corresponding intcrlami- 
nar distances are in the range 12.2-12.5 A. No other 
plateaus are seen for larger water contents. In summary, 
these results confirm the conclusions extracted from the 
preceding paragraph. 

Let us focus on the structure of the particular configu- 
ration of 27 water and 6 ethane molecules per interlam- 
inar space. This is shown in Fig. where the density 
profiles of oxygen, hydrogen, methyl, and sodium sites 
are plotted vs. the interlaminar distance from the clay 
surface. Note that oxygen, hydrogen, and sodium pro- 
files are nearly equal to those reported by Chavez-Paez 
et al. 24]. That is, we observed the same peaks at the 
same distances for the different species. The only dif- 
ference is that the water peaks are shorter due to their 
smaller content. Hence, the presence of ethane molecules 
does not practically affect the structure of water and ion 
sites. In other words, ethane molecules do not seem to 
interact very much with water molecules and ions. 

To study the orientation of ethane and water molecules 
as a function of the interlaminar position, we compute 
the statistical average of the second Legendre polynomial 
order parameter 

P 2 = (3 cos 2 (#0/2- 1/2) (8) 

for all molecules, i. In case of water, the interlaminar 
oxygen position and two angles were considered. The 
first one formed between the dipole moment vector of 
the considered water molecule and a vector normal to the 
clay layers (reference vector). The other one is formed 
between the hydrogen-hydrogen direction and the refer- 
ence vector. These are shown in Fig. |3 b) as solid and 
dashed lines, respectively. In case of ethane, the posi- 
tion of the center of mass of each molecule and the angle 
formed between the methyl-methyl direction and the ref- 
erence vector were considered. Hence, a positive value of 
Pi means a tendency for the vectors to get parallel to the 



normal of the clay sheets, whereas a negative value of Pi 
indicates a tendency for the vectors to get perpendicular 
to the reference vector. 

It is observed at the interlayer midplane that the wa- 
ter molecules show a small tendency to have their dipole 
moment vectors normal to the reference vector. Fur- 
thermore, they also show a tendency to have their H- 
H-vectors parallel to the reference vector. For water 
molecules closer to the clay surfaces, however, the H- 
H-vectors turn perpendicular to the reference vector and 
their dipole vectors become more normal to the reference 
vector. These results are in good agreement with other 
works 0,0]. There are, however, some differences of be- 
havior between the TIP4P and the MCY- water molecules 
previously studied by us ^^]. That is, for the TIP4P wa- 
ter we do not observe those few water molecules that 
escape from the midplane and stick to the siloxane sur- 
face that were observed for the MCY-water model. On 
the other hand, ethane molecules were observed to show 
a tendency to align parallel to the clay surfaces. 

To study how the confinement and ethane presence af- 
fect the ions solvation, the radial distribution functions 
and coordination numbers for sodium-oxygen sites were 
built (not shown). Here, total oxygen sites were split 
in water and clay oxygen sites so that we can distinguish 
between the different contributions. It was observed that 
all g(r) peak at 2.31 A, leading to a total coordination 
number of 5.9 and a partial coordination number for the 
water molecules of 3.7. Hence, sodium ions lose on av- 
erage more than two water molecules of their first shell 
to coordinate with clay oxygen sites. In addition, the 
total coordination number of 5.9 suggests that sodium 
ions are far from ethane molecules, indicating that ethane 
molecules are being left aside. We should also mention 
that since oxygen-oxygen separations on the clay surface 
are significantly smaller than oxygen-oxygen separations 
for waters of hydration shells around sodium ions, the 
total coordination numbers may be somewhat inflated. 
This may question the validity of some of these conclu- 
sions. 

The radial distribution functions for the methyl- 
oxygen, -sodium, -hydrogen and -methyl sites were also 
obtained and are shown in Fig.0] In case of oxygen sites, 
they were again split into the corresponding water and 
clay contributions. It is observed that all first methyl- 
oxygen peaks are situated at 3.55 A as well as the methyl- 
hydrogen peak. Such structures are reminiscent of those 
observed in hydrophobicity investigations and are a re- 
sult of water-water hydrogen bonding, since it is observed 
a lack of hydrogen atoms pointing toward the nonpolar 
solute. Note that the second methyl-oxygen peak ap- 
pears at 4.97 A, corresponding to the farthest methyl 
of the ethane molecule. This is also seen for the methyl- 
methyl distribution function, where a double peak at 4.13 
A and 5.13 A is observed. It should be pointed out that 
the peak for the methyl-sodium distribution function ap- 
pears as far as 5.67 A. This means that, in general, the 
water coordination shell that surrounds ions prevents the 



en 




FIG. 4: a) Radial distribution functions for methyl-oxygen 
sites. Dotted, dashed and solid lines correspond to water, clay 
and total oxygen sites, respectively, b) Radial distribution 
function for methyl-sodium (dotted line), methyl-hydrogen 
(dashed line) and methyl-methyl (solid line) sites. 



ethane molecules to approach them. Nevertheless, a very 
small shoulder appears at 3.25 A, meaning that a few ion 
shells were partially broken or that the number of water 
molecules were not enough to complete them. This may 
explain the 5.9 coordination number of sodium ions by 
oxygen atoms, instead of the 6.0 value we found for bulk. 

Fig. depicts some observations mentioned above, 
where a snapshot of this system is presented. Here, a 
side view shows how water molecules coordinate with 
ions in the middle of the topmost layer displacing ethane 
molecules. The same topmost layer configuration is also 
seen as a top view, where the solvation of sodium ions 
and the isolation of the ethane molecules are more evi- 
dent. This also happens in the lowest layer, although it is 
not clearly shown. Furthermore, it should be noted that 
ethane molecules are located practically in the middle of 
the interlaminar space where they are parallel to the clay 
surfaces, as previously mentioned. 

Fig. shows the density profiles of the interlaminar 
sites for the system containing 52 water and 6 ethane 
molecules. It is clearly seen that two well defined oxy- 
gen peaks indicate the formation of a double water layer. 
Hydrogen atoms show two peaks for each oxygen peak, a 
larger one close to the clay surface, and a smaller one al- 
most overlapping its corresponding oxygen peak. Sodium 
ions show three peaks, a larger one at the interlayer mid- 
plane and two smaller ones close to the clay surfaces. 
These observations are practically equal to those found 
elsewhere [Tij. Again, the only difference seems to be 
the height of the peaks, which is easily explained by the 
smaller water concentration in the interlayer. 

The orientation of water molecules is similar to the 
one layer case. It is observed a dipole moment vector 
and a H-H-vector mostly parallel to the clay surface for 
those water molecules close to it. For the water molecules 





FIG. 5: Snapshot of an equilibrated system having 27 wa- 
ter molecules and 6 ethane molecules per interlaminar space. 
Here, H are white, CH3 are light gray, O are dark gray and 
Na sites are black. Two views of the same configuration are 
shown, a side view at the top and a top view at the bottom. 
Note that for clarity, the bottom view shows only the topmost 
layer, representing the clay structure by lines. 



placed over the oxygen peak, it is obtained a positive 
second Legendre order parameter for the H-H-vector and 
even a slightly positive one for the dipole moment vector. 
This means that both tend to align perpendicular to the 
clay surface. On the other hand, the ethane molecules 
behave quite differently than in the previous case. Here, 
although they still situate in the midplane of the inter- 
laminar space, one CH3 site situates closer to the upper 
clay surface and the other one situates closer to the lower. 
Hence, and as can be seen in Fig. Elb), ethane molecules 
tend to align perpendicular to the clay surfaces. This 
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FIG. 6: a) Density profiles of oxygen, hydrogen, methyl, and 
sodium sites. They were obtained by fixing 52 water and 
6 ethane molecules, b) The corresponding second Legendre 
polynomial order parameter for water and ethane. 
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FIG. 7: Interlaminar distance as a function of the number 
of water molecules per interlaminar space. Symbols □, o, a, 
v, o and < correspond to 0, 6, 9, 12, 15 and 18 fixed ethane 
molecules, respectively. The data is obtained for 4 km of 
burial depth. The corresponding data obtained for ground 
level conditions (Fig. are presented as dotted lines. 
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FIG. 8: Interlaminar distance as a function of the number of 
ethane molecules per interlaminar space. Symbols □, o, a, v 
and <> correspond to 0, 16, 32, 48 and 62 fixed water molecules, 
respectively. The data is obtained for 4 km of burial depth. 
The corresponding data obtained for ground level conditions 
(Fig-HJ are presented as dotted lines. 
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FIG. 9: a) Density profiles of oxygen, hydrogen, methyl, and 
sodium sites. They were obtained by fixing 27 water and 
6 ethane molecules, and for 4 km of burial depth, b) The 
corresponding second Legendre polynomial order parameter 
for water and ethane. 



explains the two CH3 peaks shown in Fig. HJI a) . 

The interlaminar distance as a function of the number 
of water molecules for 4 km of burial depth is shown in 
Fig. [7| The different data series correspond, as in Fig. ^ 
to different ethane contents. This figure also shows as 
dotted lines the data obtained for ground level condi- 
tions. This is just to make easier the comparison. As 
can be seen, the same trend is observed for both condi- 
tions. The only difference is that for 4 km of burial depth 
the interlaminar distances are larger than for the ground 
level case. This difference is not very pronounced in any 
case, but seems to be larger for those systems contain- 



ing a larger number of molecules. Similarly, Fig. [S] shows 
the interlaminar distance as a function of the number of 
ethane molecules, keeping the number of water molecules 
fixed. Furthermore, and as in the previous figure, the 
data obtained for ground level is presented as dotted 
lines. Again, it is observed that the tendencies are sim- 
ilar although larger interlaminar distances are obtained 
for the largest temperature and pressure case. These 
findings are similar to those previously reported by us 

m 

Finally, the density profiles and second Legendre poly- 



nomial order parameters for the different sites are shown 
in Fig. 1^1 for the system having 27 water and 6 ethane 
molecules and for 4 km of burial depth. As can be seen, 
both, density profiles and order parameters are very sim- 
ilar to those obtained for ground level conditions (see 
Fig. Gljl. Nevertheless, the peaks are wider and shorter 
than those obtained for ground level conditions, indicat- 
ing a larger disorder of the interlayer structure. This is 
in good agrement with previous simulations |l4l 1 1 5l | and 
experiments |36). 



B. Sampling in the fj,P zz T ensemble 
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As mentioned in section III Bl this ensemble allows 
obtaining the interlaminar distance, water content and 
ethane content for the studied system at equilibrium with 
any given reservoir's conditions. For that purpose, the 
chemical potential of the species contained in the reser- 
voir at given conditions must be known. Hence, expres- 
sion (3/j,=/3fj, +hi (p/po) was used, where po is the vapor 
pressure at equilibrium with liquid water whose chemi- 
cal potential is fiQ, and p is the vapor pressure. For the 
TIP4P water model and for T=298 K and P=l atm, we 
obtained flpo = -17.4 by simulating 200 water molecules 
in a cubic box. This value is in good agreement with 
others reported for the same model [24j ■ For the studied 
burial conditions, i. e., for T=394 K and P=600 atm, we 
obtain (3/j.q = -13.4. 

For obtaining the water chemical potential, we em- 
ployed the Rosenbluth insertion method as described in 
detail elsewhere [33,|33|. For this purpose, fc o =100 oxy- 
gen trial sites are randomly chosen inside the simula- 
tion box and the quantity exp[— /3ui] is evaluated for 
each case. Here, Ui is the potential energy for a par- 
ticular oxygen site that comes only from the Lennard- 
Jones contribution. Immediately after, a favorable oxy- 
gen site is selected by assigning to each site the prob- 
ability Pi=exp[—(3ui]/W with W =J2i exp[— j3ui\. For 
this selected oxygen site, fc^=20 hydrogen configurations 
are tried while the quantity exp[— 0Uj] is calculated. In 
this case, Uj contains the contribution of the real part 
of the electrostatic potential. Hence, both hydrogen 
atoms and the charged non massive site of the TIP4P 
molecule contribute to Uj. Again, a favorable orienta- 
tion is selected by assigning to each one the probabil- 
ity Pj=exp[— f3iij]/Wh where Wh=Ylj ex P[ — P u j]- Once 
the TIP4P molecule is fully inserted, the reciprocal 
space contribution to the electrostatic potential, u r , is 
evaluated. Finally, the Rosenbluth weight factor is 
WH 2 o=W Wh exp[— f3u r ]/(k kh)- Since we perform the 
sampling from a NP ZZ T ensemble, the chemical potential 
finally reads [3i| 



(3jjL = — In 



1 



hi 



N+l 



(VW Hl 



o, 



(9) 



where (• ■ ■) is the average obtained over 100000 config- 
urations in which 10 water insertions per configuration 
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FIG. 10: Typical evolution of the interlaminar space, wa- 
ter, and ethane content. Solid lines are averages whereas 
dashed lines correspond to each interlaminar space. The es- 
tablished conditions were T=298 K, P=l atm, p/po=0.1 and 
p e =0.893P. Initial conditions were established as 11.5 A of 
interlaminar space and 10 water molecules. 



were done. 

On the other hand, we use the relationship 
P^PVid.gas+^ifiPe/®) where p° Ugas is the ideal gas 
chemical potential of the reference state, p e is the ethane 
partial pressure, and $ is the ethane fugacity coefficient 
[33j. This coefficient is 0.992 for ground level conditions 
(taken from Din j^]) and 0.536 for 4 km of burial depth 
(extrapolation of the data given by Din |37jl. 

In this way, for ground level conditions and for a reser- 
voir having p/p = 0.1 and p e = 0.893P, Fig. ITU1 shows 
how the system reaches equilibrium starting from a con- 
figuration of 10 water molecules and an interlaminar dis- 
tance of 11.5 A. For this case, it is observed that equilib- 
rium is reached after a few time steps, yielding an average 
interlaminar distance of 12.46 A, an average water con- 
tent of 27.6 molecules, and an average ethane content of 
4.88 molecules per interlaminar space. It should be men- 
tioned that other runs need a larger number of time steps 
for achieving equilibrium and so, runs must be monitored 
in order to establish the step range to average. In this 
case averages were performed in the range 5000-23000 
steps. This same procedure was then repeated systemat- 
ically for obtaining swelling curves for different reservoir 
and initial conditions. 

Fig. II II was built for ground level conditions and for a 
zero ethane partial pressure. It shows two swelling curves 
which differ on the initial conditions. One starts from an 
interlaminar distance of 11.5 A and 10 randomly placed 
water molecules. The other one starts from an interlam- 
inar distance of 16.0 A and 60 water molecules. It can 
be seen that for p/po < 0.2 both curves coincide, yield- 
ing interlaminar distances and water contents lower than 
12.5 A and 42 molecules, respectively. In particular, for 
p/po = 0, no water molecules and 10.22 A of interlam- 
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FIG. 11: Fig. 1111 Interlaminar distance and number of water 
molecules per interlaminar space as a function of water vapor 
pressure and for ground level conditions. The ethane partial 
pressure was set to zero. Symbols □ and o correspond to 
the data obtained by simulations and starting from 10 and 
60 water molecules, respectively. Symbol ■ corresponds to 
experimental data reported by Brindley and Brown |3ql . 



inar distance were obtained. For p/po > 2, the curves 
separate producing two different plateaus. One corre- 
sponds to the formation of a single water layer and the 
other one to the formation of a double layer. This first 
plateau yields values of the interlaminar distance in the 
range 12.34-12.73 A, having 34.6-47.5 molecules of wa- 
ter content. The double layer plateau yields 14.82-15.57 
A and 68.0-82.1 molecules, respectively. The interlam- 
inar distance values are in good agreement with others 
reported by experiments and simulations 0, 0, 0, |35( . 
The water contents agree with the ones reported for the 
TIP4P water model and obtained by /iVT simulations 
[24j . Nevertheless, they are larger than others [34[. 

The presented swelling curves are similar to those re- 
ported by Hensen and Smit [34|. Both predict a single 
layer for p/po < 0.2 and two layers for the range p/po = 
0.3-0.55. For larger p/po, our data still predict the two 
configurations to be stable and theirs show that the sin- 
gle layer is not stable anymore. This may be the main 
difference between both results. On the other hand, for 
a similar water-clay model and sampling from a \xVT en- 
semble, Chavez-Paez et al. |2J| also found that both con- 
figurations are stable, in good agreement with our data. 
Finally, swelling experimental data taken from Brindley 



FIG. 12: Interlaminar distance, number of water molecules, 
and number of ethane molecules per interlaminar space as a 
function of water vapor pressure. The ethane partial pres- 
sure, p e , was set to 0.893P and ground level conditions were 
established. Dotted lines correspond to p e — 0, which were 
included to make easy the comparison. 



and Brown's work |35| were also included in the inter- 
laminar distance graph, in order to be compared with 
our results. As can be seen, the agreement is very good. 

A similar study was made to determine the influence of 
a ethane rich reservoir on the interlaminar space. For this 
purpose, a partial ethane pressure of 0.893P was set for 
the same ground level conditions. The obtained results 
are shown in Fig. ^| together with the curves already 
shown in Fig. ^] It is seen that, in general, they are 
very similar, i. e. there is not a very pronounced inter- 
action between the ethane atmosphere in contact with 
the clay-water system. In fact, for p/po > 0.4 the largest 
ethane content observed was 0.12 molecules per interlam- 
inar space. This is equivalent to just one ethane per 350 
water molecules. Furthermore, for the double-layer cases 
the ethane content is even lower, and this is why we did 
not include it in the plot. Moreover, the ethane propor- 
tion decreases with increasing the vapor pressure yielding 
zero for p/po > 0.8. Consequently, those points obtained 
for p/po > 0.4 practically show the same interlaminar 
distance and water content previously found. 

Nevertheless, for p/po < 0.4 a clear contribution of 
ethane to the interlaminar space can be seen. This con- 
tribution turns more important by decreasing the vapor 
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FIG. 13: Interlaminar distance, number of water molecules, 
and number of ethane molecules per interlaminar space as 
a function of water vapor pressure and for 4 km of burial 
depth. Symbols □, o and a correspond to p e — 0.0P, 0.107P 
and 0.268P, respectively. 



pressure. This means that some water molecules are 
being replaced by ethane molecules, keeping the inter- 
laminar distance almost constant. The exchange ratio 
is close to 3:1, as found in the preceding section. For 
the extreme case of p/po=0, the whole water layer is re- 
placed by a layer of 16 nearly parallel ethane molecules, 
which leads to a practically equal interlaminar distance. 
It should be mentioned that, for this particular case, the 
simulation was started from 12.5 A of interlaminar space 
and 10 ethane molecules, instead of 11.5 A and 10 water 
molecules. If started from this last initial condition, the 
system quickly loses its water content producing a sta- 
ble dehydrated state having zero ethane molecules. We 
also observed that this stable dehydrated state is destabi- 
lized by any small but enough vapor pressure. This indi- 
cates that the presence of a small amount of water, which 
widens the interlaminar space and solvate ions, aids the 
entrance of ethane molecules. This finding agrees with 
experimental evidence reported by Barrer |38| . 

Similarly to Fig. ^1 Fig. 1131 shows swelling curves but 
for 4 km of burial depth. Here, three different ethane 
partial pressures were considered. These are p e = 0.0P, 
0.107P and 0.268P. For p e — 0, some differences with re- 
spect to the ground level case can be observed. 

On the one hand, for p/po—1.0 neither a single layer 



nor a double layer configuration are stable. For this 
particular case, no matter what the established initial 
conditions, the interlaminar space and water content 
monotonously increase as the simulation evolves. They 
were ended after reaching 180 water molecules and, there- 
fore, we assume that a full hydration of the clay sys- 
tem was obtained. A similar behavior was observed for 
p/po=0.8 and for an initial configuration of 60 water 
molecules. When starting this last system from an initial 
configuration of 10 water molecules, however, the system 
reaches equilibrium yielding a single water layer. This is 
in good agreement with the results reported by Chavez- 
Paez et al. 0], where just a stable crowded single water 
layer configuration was obtained for burial conditions. 
Nevertheless, they obtain this behavior for p/po=1.0, in- 
stead of p/pa=0.8. Their employed conditions were 353 
K and 625 bar, which correspond to depths of 4.16 and 2 
km, in terms of lithostatic and thermal gradients, respec- 
tively. Hence, it seems that the difference of temperature 
employed in both cases may explain the different values 
oip/po that produces a stable single layer and an instable 
double layer. 

On the other hand, it is observed that both curves 
predict the same single water layer for p/po<0.3. Con- 
sequently, the double layer configuration is stable just 
in the range p/po — 0.4-0.6. This result agrees with 
the experimental finding of a dehydration of the inter- 
laminar space from a double water layer to a single one 
as the burial depth increases [3(|. In addition, it also 
agrees with previous simulations of a similar system but 
for the MCY water model, where the stability of the two 
layer configuration was found to decrease with increasing 
burial depth |l5j . 

Additionally, the interlaminar distances for single and 
double layer configurations are similar to those corre- 
sponding to ground level conditions. The water con- 
tent, however, is slightly lower, suggesting that water 
molecules occupy a larger effective volume. This finding 
also agrees with the results reported for the MCY water 
model [3. 

For p e = 0.107P and 0.268P, a similar effect to the one 
found for ground level conditions is seen. That is, some 
water molecules are replaced by ethane molecules, and 
this exchange is enhanced for lower water vapor pres- 
sures. The extent of this exchange, however, is much 
larger under burial conditions. This is mainly explained 
by the higher ethane pressure, which leads to a larger 
ethane activity. Naturally, the water-ethane exchange is 
larger for p e — 0.268P than for 0.107P for all cases, as 
expected. It is observed that even for p/po =0.6, 0.4 and 
0.9 ethane molecules enter the interlaminar space for the 
single layer case and for p e = 0.107P and 0.268P, respec- 
tively. This is approximately equivalent to one ethane 
molecule per 100 and 45 water molecules. For the dou- 
ble layer cases, the water-ethane exchange is not so pro- 
nounced. For p e =0.268P and p/po =0.6, only one ethane 
molecule per 220 water molecules were detected by aver- 
aging a large number of /j,P zz T ensemble configurations. 
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In spite of the important exchange of water by ethane 
molecules found for p/j>o<0.6, several water molecules 
still remain in the system even for p/po as low as 0.05. 
For this case and for p e =0.268P, an average of 20.8 water 
molecules were found to be strongly attached to cations, 
impeding the entrance of more than 7.8 ethane molecules. 
As for the ground level case, when setting p/p —0, wa- 
ter is forced to leave the system and a single layer of 16 
ethane molecules is obtained. Again, this layer is not ob- 
tained if the simulation is started from an initial config- 
uration close to the dehydrated state. This means that 
a dehydrated clay containing just ions in its interlaycr 
space and in contact with a dry ethane rich reservoir 
does not allow the entrance of ethane molecules. 



IV. CONCLUSIONS 

Hybrid Monte Carlo simulations in the NP ZZ T and 
/j,P zz T ensembles were used for studying the Na- 
montmorillonitc hydrates under ethane rich reservoirs. 
From the NP ZZ T study, interlaminar distances as a func- 
tion of the water and ethane content were obtained. 
These curves show clear plateaus for lower ethane con- 
tents and mainly for water contents consistent with the 
formation of a single water layer. 

On the other hand, from the NP ZZ T ensemble we also 
studied the interlayer structure. Here it was observed, 
for systems containing few ethane molecules and water 
enough to complete a single layer, that ethane mainly 
situates close to the interlayer midplane and adopts a 



nearly parallel arrangement to the clay surface. Further- 
more, it was observed that the water-ion-clay interactions 
are practically unchanged by the presence of the hydro- 
carbon, indicating that this component is being left aside. 

In addition, the \xP zz T ensemble allows us determin- 
ing the interlaminar distance and water-ethane content 
for any specific reservoir's conditions. Hence, several 
swelling curves as a function of water pressure were ob- 
tained. For ground level conditions, they indicate that 
ethane molecules do not practically enter the interlami- 
nar space, unless one establishes a very low water pres- 
sure. If this is the case, we also observed that water 
molecules are flushed out by this entrance, leading to an 
almost constant interlaminar space close to 12.5 A. 

For the studied burial depth conditions, it was ob- 
served that the exchange is enhanced, and that the 
ethane molecules are capable of entering even for rela- 
tively large water pressures. On the other hand, this ex- 
change is not very pronounced for a double layer config- 
uration. Furthermore, it was observed that for very low 
water pressures, the system holds water enough to keep 
the ions' shells practically intact. Finally, our results also 
indicate that a dehydrated and collapsed clay system is 
not capable to swell and absorb ethane molecules if no 
water vapor is present in the atmosphere. 
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